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Abstract 

I present and analyze a quadratically convergent algorithm for computing 
the infinite product FI^L 1 (f — tx n ) for arbitrary complex t and x satisfying 
\x\ < 1, based on the identity 

TJ(l-tx n ) = V ( } - 

J-J- 1 ' ^ (1 - x)(l - x 2 ) ■ ■ ■ (1 - x m ) 

n=l m=0 v Jy ' v ; 

due to Euler. The efficiency of the algorithm deteriorates as |x| \ 1, but much 
more slowly than in previous algorithms. The key lemma is a two-sided bound 
on the Dedekind eta function at pure imaginary argument, r](iy), that is sharp 
at the two endpoints y = 0,oo and is accurate to within 9.1% over the entire 
interval < y < oo. 



2000 MATHEMATICS SUBJECT CLASSIFICATION: Primary 33F05; Sec- 
ondary 05A30, 11F20, 11P82, 33D99, 65D20, 82B23. 

KEY WORDS: Euler's partition product, g-series, g-product, Dedekind eta func- 
tion, numerical algorithm. 



1 Introduction 

The function 



R(t,x) = Y[{l-tx n ) , 



n=l 



defined for complex t and x satisfying |x| < 1, was first studied by Euler |13| and 
has numerous applications in combinatorics, number theory, analytic-function theory 
and statistical mechanics. The case t = 1 is equivalent to the Dedekind eta function 



7]{t) 



^mr/12 



R{l,e 



2nir\ 



X2) 



which is a modular form ]5], ^TJ and plays a central role in the enumeration of partitions 
||. H ^TJ and sums of squares | 2lf . The case t = — 1 is related to t = 1 via the trivial 
identity 

Both of these cases are related to theta functions 0, ID, 11, via the identities 

oo oo 

R(l,x) = Y[{l-x n ) = {-l) m x m{3m+1)/2 (1.4) 

n=l 

oo 

R(i,xf = n(i 

R(l,x 2 ) 2 



n=l 
oo 



m=— oo 



X 



n\3 



1 — X 



2n 



m=0 

oo 



n=l 



— X 



2n-l 



R{1,X 



2\5 



i?(l, x) 2 i?(l, X 4 ) 2 



R(l,x) 

oo 

= JJ(l-x 2n )(l + a; 2n - 1 



x m(m+l)/2 



m=0 



n=l 



R(l,x) 
R(l,x 2 



n 

71=1 



1 - x r > 



X" 



^(-l) m (2m + l)x m(m+1)/2 (1.5) 

(1.6) 

l= — oo 
oo 

£ (-l) m x m2 (1.8) 



due to Euler, Jacobi and Gauss, which have spawned a plethora of modern extensions 
fl6| , [261 , H; HH) S IH> HI- Additional cases of the function R(t,x) arise in the 
celebrated Rogers-Ramanujan identities P, ffl 



n=0 v 



1 



-X 5n+1 )(l - X 5n+i ) 



^-i (1 - aO(l 



;r 



m=0 



nn 



(1 - z)(l - X 2 ) • ■ ■ (1 -x m ) 

^,m(m+l) 



71=0 



(l_ x 5n+2)( 1 _ a .5n+3) 



^ (l-:r)(l 



m=0 



(1 - x)(l - X 2 ) ■ ■ ■ (1 -x m ) 



(1.9) 
1.10) 
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which have numerous combinatorial consequences [0] and which play a key role in 
Baxter's solution of the hard-hexagon problem in statistical mechanics || [7]]. (See 
also |TJ], 0] for many related identities.) Finally — and this was the initial 



motivation for the current work — the cases t = ±1 and t = ±to, where u is a 
cube root of unity, arise in Baxter's solution for the chromatic polynomials of large 
triangular lattices || || . To determine the limiting curves of chromatic roots for these 
lattices, it is necessary to compute R(t,x) to high precision for complex x, including 
points x very near the unit circle |Tj| . 

The numerical computation of R(t,x) clearly becomes delicate when \x\ | 1. 
Surprisingly, there seem to be very few treatments of this problem in the litera- 
ture |34|, |35|, [15|, |IJ, and the algorithms employed there are only linearly convergent; 
moreover, these authors (with the exception of Gatteschi |TJ|) considered almost ex- 
clusively the case of real t and x. My purpose here is to propose and analyze a 
quadratically convergent algorithm for computing R(t, x) for arbitrary complex t and 
x satisfying \x\ < 1, based on the identity 

00 00 (f\m m(m+l)/2 
R(t, X ) = T7(l-tx n ) = Y~ 9 , r (1.11) 

due to Euler.[] In the course of this analysis, I will obtain (Corollary |2.4| ) a two-sided 
bound on R(l,x) for < x < 1 (and thus on r](iy) for < y < oo) that is sharp at 
the two endpoints x — 0, 1 and is accurate to within 9.1% over the entire interval; 
this bound is perhaps of some modest independent interest. 

Of course, for the special case t = 1 one may employ an even faster algorithm 
based on using the modular transformation law for the Dedekind eta function 0, 
[5], [□], 31] to move x away from the unit circle, followed by evaluation of the 



quadratically convergent sum ( p..4|) f| Moreover, the case t = —1 can be reduced to 
t = 1 via ( |1.3D . But for t 7^ ±1 no such identities are known. 

The plan of this paper is as follows: In Section ^] I formulate and prove the 
properties of the function R(t, x) that will be needed in the sequel. In Section |3| I 
obtain bounds (both a priori and a posteriori) on the rate of convergence of the 
algorithm defined by ( |1 . 1 1| ) . Finally, in Section |], I briefly compare this algorithm to 
other algorithms that have been proposed [Tj| [I[ for computing R(t, x). 



1 For a proof of ( pTTl| ), see e.g. §, p. 19, Corollary 2.2], @ p. 34, Lemma 4(a)] or || PP- 22-23]. 

2 Surprisingly, I have been unable to find in the literature any discussion of such an algorithm. The 
details of its implementation — in particular, how to find an appropriate modular transformation 
— may not be entirely trivial. 
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2 Properties of R(t, x) 



We shall assume throughout this note that |x| < 1, even if it is not explicitly 
stated. Likewise, when we write x = e~ 7 , we shall assume that Re 7 > 0. 

2.1 Elementary properties 

We begin by noting some elementary properties of R(t,x): 

1) R(t, x) is a jointly analytic function of t and x for (t, x) G C x O, where D is 
the open unit disc. For fixed R(t, x) is an entire function of t of order 0, with 
simple zeros (for x 7^ 0) at t = x~ l , x~ 2 , .... 

2) By splitting off the first term in the product ( |1 . 1|) , we obtain the functional 
equation 

R(t,x) = (1 - tx) R(tx, x) , (2.1) 

from which Euler's formula ( |1 . 1 1| ) can easily be derived by comparing coefficients of 
powers of t. 

3) Let wbea primitive mth root of unity, and use the identity nj^o^l — u " ' z ) = 
1 — z m ; we find 

m— 1 

Y[ R{uH, x) = R{t m ,x m ) . (2.2) 
3=0 

The special case m = 2, t = 1 is (|L~3"1). 

4) By splitting the product (|1 . 1| ) according to residue classes modulo m, we obtain 



R(t,x) = Y[ R(tx j - m , x m ) . (2.3) 
i=i 

This formula permits the determination of the asymptotic behavior of R(t, x) as x 
approaches an mth root of unity, once the asymptotic behavior as x — > 1 is known. 

5) We have the trivial upper bound 

\R(t,x)\ < R(-T,\x\) (2.4) 

whenever \t\ < T. 

6) We have the trivial lower bound 

\R(t,x)\ > R(T,\x\) (2.5) 

whenever \t\ < T < (The condition T < is of course best possible, since 

R(t,x) vanishes at t = x" 1 .) 

7) Finally (and most importantly), let us take the logarithm (principal branch) of 



the defining equation (|1 . 1|) , expand log(l — tx n ) in Taylor series, and interchange the 
absolutely convergent summations; this yields the useful representation as a Lambert 
series 

]ogR(t,x) = -Etj-T^' ( 2 - 6 ) 

k=l 

valid whenever \x\ < 1 and \tx\ < 1. We shall use this representation repeatedly. 
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2.2 Elementary bounds 

Bounding the denominator of ( |2.6| ) using |1 — x k \ > 1 — > 1 — \x\, we obtain: 

Lemma 2.1 Whenever \x\ < 1 and \tx\ < 1, we have 

\\og R(t,x)\ < f u (2.7) 

(where the principal branch of the logarithm is taken) and hence 

(1 - |tx|) 1/(1 - |xi) < \R(t,x)\ < (1 - . (2.8) 

This is a crude bound that does not exhibit the correct behavior as |x| — > 1, but we 
shall use it starting point for further refinements. 



First we need a slight extension of Lemma |2.1| for the special case t = 1. For 
< x < 1, define 

00 i k 

S(x) = -logics) = EfcT3^' ( 2 - 9 ) 



fc=i 



so that 

S'lx) 



fc=i 



s „ (i) = E (*-i)^ + i)^ ( , n) 

k=l ^ > 

By using 1 > 1 — x h > 1 — x in the denominator, we obtain the trivial bounds: 
Lemma 2.2 For < x < 1, we have 

-log(l-x) < S(x) < < ^—^ (2.12) 

1 - t^~i - £ jdw (2 - 13) 
3 < i + < g-fai < ^+ . 2 :/ 2 ... < 3 



(1-x) 2 (1-x 2 ) 2 - v y - (1-x) 5 (1 -x) 3 (l -x 2 ) 2 - (1-x) 5 

(2.14) 



2.3 Case t = 1 



Now we improve these bounds by using a deep fact: the transformation properties 
of the Dedekind eta function under the modular group || ||, [□], ^T], [311 . All we 



5 



need, in fact, is a special case of the modular transformation law, namely the one for 
inversion r — > — 1/r: 



i?(l,e"^) = ^ Vi exp^- - — j i2(l,e-*/*) (2.15) 

for Kez > 0.0 This allows us to control the behavior near x = 1 (z — > 0) in terms of 
the (trivial) behavior near x = (z — > +oo).[] Indeed, from ( |2.15 ) and the regularity 



of R(l,x) near x = 0, one immediately deduces the sharp asymptotic formula 



tt 2 1 1 , . 7 



logit!(l,e^) = -— - - log 7 + §log(27r) + ^ + 0(e^ h) (2.16) 

as 7 — > 0; moreover, an explicit quantitative bound on the 0(e~ 47r2//7 ) term can easily 
be extracted from Lemma |2.1| . 

For later applications we need also a quantitative error bound valid for real 7 in 
the entire interval < 7 < 00. Let us define 

f(z) = logflCL.e-*") + ^ - ilog^l + i^ , (2.17) 

so that 

/(1/z) = lo gj R(l,e- 2 ^) + ^ - ^log(l + ^ 2 ) . (2.18) 

Then the transformation law fl2.15|) tells us immediately that f(z) = f(l/z), and 
indeed we have: 

Proposition 2.3 For < z < 00, we have: 

(a) f(z) = f(l/z) 

(b) ]im f(z) = and lim f(z) = 

2J.O z—*+oo 

(c) < f(z) < f(l) = f - J log 2 + log 77^) « 0.0866399 

(d) f (0) = tt/12, /'(«) > /or < 2 < 1, /'(l) = 0, and f(z) < /or z > 1 

(ej /"(0) = —1/2, and taere exists z* > 1 stzca taat /"(z) < /or < z < z* and 

f"M = 0. 



3 There are a number of proofs of (2.15). The simplest uses the Poisson summation formula 



applied to Euler's pentagonal number theorem ( |l.4l) 21 Section 3.3]. Another proof, due to Siegel 



uses the Cauchy integral formula J 5[ Section 3.2] fl l|, Section VIII. 3]. Proofs of the full modular 
transformation law are given in Ji, Sections 3.3-3.6 and pp. 190-195], Q Sections 3.1-3.3 and 
4.1-4.2], f| Chapter 9], and §, pp. 82-85]. 

4 Using the full modular transformation law, one can control in an analogous way the behavior 
of R(l,x) near any point x = e 2mh ^ k (h, k G Z) of the unit circle: see e.g. j| Chapter 5]. 
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Proof. We have already proven that f(z) = f(l/z), so we can use ( p,17| ) and ( f2.18| 



interchangeably as formulae for f(z). The limiting values of / and its derivatives at 
z = can be read off fl2.18| ). 



To prove f(z) > for < z < oo, it suffices to prove it for < z < 1. Using 
( |2.18| ), we make the following crude bounds: 

f(z) = -S(e~ 2 ^) + ^ - ^log(l + z 2 ) (2.19a) 

e- 2n / z ttz 1 

£ -713^7^ + 12 "1 108(1 + 2 > (2 ' 19b) 



e 



7T2; 


1 


12 


4 




TTZ 


- + 




12 " 






- + 


(- 


\ 12 



> ~ + ( - - 4 ) * (2-19c) 



where we have used ( |2.12| ) and the fact that < z < 1. So we need only show that 

1 <fe4V-^ (2.20) 



(1 - x) 2 \12 A J \ logx 

for < x < e~ 27T . But — xlogx/ (1 — x) 2 is an increasing function of x for < x < 1, 
and its value at x = e" 2 *" is 27r e - 2 7(l - e^) 2 « 0.011777 < (vr/12 - 1/4)(2tt) « 
0.074138. 

Next let us prove that there exists e > such that f"(z) < for < z < 1 + e. 



Differentiating (2.18j) twice with respect to 2, we obtain 



/"(a) = + e - 2 -A (S "( e -2 7 rA) _ l^ e - 4 ^5"(e- 2 ^) ' 



^ 3 y v v ' 2(1 + ^ 2 ) 2 ' 

(2.21) 

From (|2.13| ) / (|2.14|) we have S'(x) > 1 and S"(x) > 3, so the first two terms in ( 2.21) 



are < for < z < 7T, and the third term is < for < z < 1. This proves the 
claim. 

We have just proven that f'(z) is a strictly decreasing function of z on < z < 1+e. 
From f(z) = f(l/z) it follows that f'(l) = 0. Therefore f'(z) > for < z < 1; by 
/(z) = /(l/z) it follows that f (z) < for z > 1; and thus /(z) < /(l) for all z. 

Finally, it is not possible that f"(z) < for all z, as this would imply that f(z) < 
for some z G (1, oo). So we can define z* > 1 to be the smallest z such that f"(z) = 0. 
■ 

Remark. Numerical calculations show that /" has a unique zero, which is located 
at z* ~ 1.974174. But we shall not bother to prove this. Graphs of f(z) versus z and 
log z are shown in Figure the latter shows the z 1/z symmetry more clearly. 



47T 2 ^ 



Proposition |2.3| can be rephrased by defining 

i2o(l,a;) = e^'^°^ (l + -^—r ) , (2.22) 



v (logx) 2 

which we interpret as an "approximate" version of R(l,x). We then have: 
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Corollary 2.4 For < x < 1, 

e ^/(6io g x) < ^(x^) < e{1,x) < CRq(1,x) (2.23) 
where C = e /(1) = e^ -1 / 4 ^) « 1.090504. 

In other words, we have a two-sided bound on R(l, x), in which the lower bound is 
sharp at the two endpoints x — 0, 1 and is accurate to within 9.1% over the entire 
interval < x < 1. We shall frequently use the lower bound of Corollary [T4] in the 
form 



9, / 4vr 2 \ 1 
#(l,e~ 7 ) > e- 7r/67 (l + — J > e-"/ 67 (2.24) 



for 7 > 0. 



2.4 Case t = -1 



We can now handle the case t = — 1 by using (|1.3|) to relate it to t = 1. From 



(" p. 161 ) and (|1.3|) we obtain the sharp asymptotic formula 



logi?(-l,e^) = ^ - Il 0g 2 + ^ + 0(e-^) (2.25) 

as 7 — > 0, where again a quantitative bound on the 0(e _7I " 2//7 ) term can easily be 
extracted from Lemma I2TTI. Moreover, we can obtain a quantitative error bound valid 



for real z in the entire interval < z < oo. Let us define 

g(z) = f(V2z) - f(z/V2) (2.26a) 

= l„ gfl (-l, e -^) - + ilog(^±|) ■ (2.26b, 



It follows immediately from Proposition |2^ that 
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-0.03 1 -0.03 1 

Figure 2: Graphs of g(z) versus z and log z. 



Proposition 2.5 For < z < oo, we have: 

(a) g(z) = -g(l/z) 

(b) \img(z) = and Km g(z) = 

zjO z^+oo 

(c) g\z) < for l/y/2 < z < ^2 

(d) g(z) > for l/y/2 < z < 1, g(l) = 0, and g{z) > for 1 < z < V2 

(e) \g(z)\ < /(l) w 0.0866399 for < z < oo 

Remark. Numerical calculations show that g' vanishes when (and only when) 
±log2 ~ 1.180158, i.e. z or 1/z m 3.254889, and that the maximum value of \g(z) \ is 
w 0.0251707. It follows that R(-l,x) differs from 

(7T 2 \ 1 / 4 
1+ P°ff ) (2.27) 

by less than 2.6% over the entire interval < x < 1. Graphs of versus z and 
log 2 are shown in Figure 

2.5 Asymptotics of R(t, x) for General t 

Finally, let us discuss briefly the asymptotics of R(t, x) as x — > 1 when t is fixed 
with \t\ < 1 (or more generally varies within a compact subset of the open unit disc). 
Let us write x = e~ 7 with Re 7 > and study the behavior as 7 — > 0, using the 
representation (|2.6|) . We have 

k 1 00 ry 

E^Tr- 1 (2.28) 



1 — x k e fc7 — 1 ' ml 

m=0 
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where B m is the rath Bernoulli number; this series is absolutely convergent for |7| < 
27r/fc.[] Inserting this into (|2.6|) and formally interchanging the order of summation, 



we obtaii 



where 



log R(t, 



00 R 

— r Li 2 - 
ra! 



k=l 



(2.29) 



(2.30) 



is the polylogarithm function f24|. However, because the radius of convergence of 
( j2.28| ) is nonuniform in k and tends to zero as k —>■ 00, it is reasonable to expect that 



the series ( p.29|) is not convergent but is only asymptotic. One further expects that 
this asymptotic expansion should hold uniformly as t varies within a compact subset 
of the open unit disc. All these expectations are true [36]. What is perhaps more 
surprising is that the expansion ( |2.29| ) holds also for t on the unit circle, except at 
the point t — 1. Indeed, under suitable restrictions on arg7 it holds in a much larger 
domain of the complex t-plane, which in the most favorable case (7 real and positive) 
encompasses the entire complex t-plane except for a cut along [l,oo). These results 
will be reported elsewhere ||36|| . For real 7 > and < t < 1, the expansion ( ]2.29|) 
was proven some years ago by Moak [28], Theorem 3].[] For real 7 > and t < 1, 
the expansion ( p. 29 ) and some generalizations thereof have recently been proven by 
Mcintosh P7| . For real 7 > and t £ C\ [l,oo), the expansion ( ]2.29| ) has been 
proven by Prellberg H30| , Lemma 3.2]. All these works use the Euler-Maclaurin sum 
formula. Our approach [R6J, by contrast, uses complex integration. 

For real 7 > and < t < 1, we can use the method just sketched to obtain a 
two-sided bound on R(t, e -7 ) that incorporates the first two terms of the expansion 
( p.29| ). For z > we have the elementary inequalities^ 



1 

z 



1 

2 £ e ~ 



1 1 , 



1 



< 



r z/2 l 

< - 

Z Z 



(2.31) 



Setting z = k-f and inserting these bounds into (|2.6|) , we obtain: 

Proposition 2.6 For < t < 1 and 7 > 0, we have 

- log R{t,e~ r ) < 7 " 1 Li 2 (te~ 7/2 ) (2.32a) 
< 7~ 1 Li 2 (t) (2.32b) 

5 See e.g. 0, equation (6.81)]. 

6 See also jl2], p. 58, exercise 2] and [0 Theorem 4] for this formula. 

7 Equation (4.2) of ^8| contains a misprint: there should be a minus sign before the integral. 
Correspondingly, in equation (4.3), the minus sign before the integral should be a plus sign. 

8 The first two inequalities can be derived from tanh(z/2) < z/2; the third can be derived from 
z/2 < sinh(z/2); and the fourth is trivial. Note that all of these bounds, except the last, capture 
the first two terms of the Laurent series for l/(e 2 — 1) around z = 0. 
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and 



\ogR{t,e-" / ) > 7 - 1 Li 2 (te" 7 ) - |log(l-te- 7 ) 
> 7~ 1 Li 2 (t) + ilog(l-t) 

It is worth remarking that, even for t — 1, the bounds ( 2.32a| ) and ( |2.33a ) capture 



(2.33a) 
(2.33b) 



the first two terms of the asymptotic expansion ( [2.16| ), i.e. they get the correct log 7 
term. 



One application of Proposition \LQ is to bounding the partial product 

R(l,x) 



k=l 



— X 



R(x n ,x) 



(2.34) 



when < x = e 7 <1 (and we will usually take n to be of order 1/7). Inserting the 
lower bound ( |2.24| ) on R(l,x) and the upper bound ( |2.33b| ) on R(x n ,x), we obtain: 



Corollary 2.7 Let 7 > 0. Then 



lid- 



-k-y\ 



> exp 



k=l 



Li 2 (e-™ 7 ) - tt 2 /6 

7 



:i - e -^) 1/2 



(2.35) 



In particular, for n < (log2)/7 we have 

-(lo g 2)72-^/12 



k=i 



— e 



> exp 



7 



'1 - e-" 7 ) 1 / 2 



(2.36) 



Here ( |2.36| ) follows from ( |2.35| ) and the well-known fact [23, |25 

7T 2 (log2) 2 



Li 2 (l/2) 
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(2.37) 



Mcintosh I23] has recently obtained a complete asymptotic expansion of the partial 
product rifc=i(l — te~ fc7 ) for n = /i/7 (/z fixed, real 7 j 0) and either t = 1 or t < 1. 

3 Numerical Computation of R(t, x) 



In this section we discuss the use of Euler's formula 

_^n a ,n(n+l)/2 



R(t,x) = V- ^~ 

K ' ; ^ 1 -x 1 



n=0 



(1 - x)(l -X 2 ) ■■■(!- X n ) 



(3.1) 



to compute R(t,x) for complex t and x satisfying |x| < 1. We shall give two types of 
bounds on the error committed by truncating the series ( |3.1| ): 



(a) an a priori bound in terms of \t\ and \x\ alone; and 
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(b) an a posteriori bound, based on the actual computed value of the last included 
term. 

We shall also give some guidance about the needed numerical precision in intermediate 
stages of the calculation, by comparing the largest term in the sum to the final answer. 
We use the following definitions: 

f—f\n x n(n+i)/2 

• The nth term: a n = - — — — 

(1 -x)(l -x 2 )---(l 



— x r 

N-l 

• The partial sum after N — 1 terms: Sn = ^ a n 

n=0 

oo 

• The remainder after N — 1 terms: Rn = XI 

n=N 

• The absolute error after N — 1 terms: Ajy = \Rn\ 

• The relative error after N — 1 terms: 5^ = \R^/R(t,x)\ 

• The modified relative error after N — 1 terms: 5' N = \Rn/Sn\ 

Clearly S' N /(1 + 5' N ) < 6jy < S' N /(1 — 5' N ), so the two types of relative error are 
essentially indistinguishable when 5n,5' n <C 1. 

Lemma 3.1 If \x\ < 1 and \t\ \x\ N+1 < 1, then 

_°° \f\N |„|iV(JV+l)/2 

y | rx n(n+l)/2| < El IgL . (3.2) 

^ 1 1 ~ 1 - t \x\ N+1 v ; 

n=N I M I 

Proof. Bound the sum by a geometric series, using 

: |t||x n+1 | < \t\\x\ N+1 (3.3) 



n=N 



^n+l x (n+l)(n+2)/2 



-j-Tij,n(n+l)/2 



for n > N. 



Lemma 3.2 If \x\ < e 7 (j > 0), then 



11(1 -«*) 



k=l 



> JJ(l-\x\ k ) > Y[(l-\x\ k ) = R(l,\x\) > e"" 2/67 . (3.4) 



k=l k=l 



Proof. An immediate consequence of Corollary |2 

Remark. An improved bound on the partial product n]? =1 (l — x k ) can be obtained 
from Corollary |2.7| ; it is advantageous when 717 ^> 1. 
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Proposition 3.3 Suppose that \x\ < e 7 wif/i 7 > 0. 

UlN e 7r 2 /67-AT(Af+l)7/2 



faj < e^ 1 ^, then A 



iV 



n=Af 



< 



u, e _ (iV +l) 7 



n=N 



N 



(c) If\t\ < e 1 , then 8 



\R(t,x) 



E a n 

n=N 



N 



\R(t,x) 



< 



< 



3 7T 2 /37-iV(jV+l)7/2 
1 _ e -(iV+l) 7 ■ 



3 7r 2 /3 7 -Af(Af+l)7/2 



1 



1 - e -( N+1 ^ 1 - lile" 7 ' 



Proof, (a) is an immediate consequence of Lemmas |3.1| and |37|. (b) follows from (a 
together with the bound \R(t,x)\ > R(l, \x\) > e _7r Z 67 from ( |2.5| ) and Corollary 
(c) follows from (b) and (|2.1| ). ■ 



Corollary 3.4 Let K > 0, and suppose that \t\ < 1 and \x\ < e 7 (j > 0). 



71 



2K 



(a) If N > \ — H , then A N < e 



-K 



3 7 



7 



(b) IfN > 



'2tt 2 
37 1 



2K 

7 



then 5n < e 



-K 



Proof. Since K > 0, we have A^7 > and hence 

e -Ny/2 e -ir/2V3 



-(Af+1)7 



< 



0.482426 < 1 . 



(3.5) 



Now vr 2 /67 - 7Y 2 7/2 < -K in case (a), and tt 2 /37 - N 2 ^/2 < -K in case (b). The 
result then follows from Proposition |3.3| (a,b) . ■ 

Please note that the bound in Proposition |3.3| (a) is asymptotically within 9.1% 
of being sharp when < x = e~ 7 < 1 and N ^> I/7 (and in this case is moreover 
asymptotically sharp as 7 j 0); but it is overly pessimistic in other cases, because the 
denominator (1 — x)(l — x 2 ) • - • (1 — x n ) is not really as small as Lemma |3.2| says it 
could be. Likewise, the bound in Proposition |3l|(b) is asymptotically (almost-) sharp 
when, in addition to the above conditions, we have t — 1; but it is overly pessimistic in 
other cases, because \R(t,x) \ is not really as small as the bound \R(t, x)\ > R(l, \x\) 
says it could be. 

It is thus of some value to provide an a posteriori bound on the truncation error 
that is more realistic, when x ^ (0, 1), than the a priori bound; such a bound can 
be used a stopping criterion in the numerical algorithm. We need the following 
elementary observation: 
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Lemma 3.5 If \x\ < e 7 with 7 > 0, then 



\t\ \x\ n 
11 - x n \ 



< 



\t\ \x\ 



\x\ 



< 



\t\ e 



-rvy 



1 - e" n T 



(3.6) 



Lemma |3.5| tells us that, at least for < x < 1, the terms a n increase in magnitude 
until \x\ n ~ 1/(1 + i.e. n ps [log(l + |£|)]/7, and then decrease. (For general 
complex x, the terms will sometimes increase up to this point, i.e. for those n for which 
argi" ~ mod 2n. How often this occurs depends on the Diophantine properties of 
argx.) We can use Lemma [3]5] to bound the tail of the sum by a geometric series: 

Proposition 3.6 Suppose that \x\ < e™ 7 (7 > 0) and N > [log(l + |t|)]/7- Then: 

\t\e- N ^ 



(a) A 



N 



(b) 5' T 



n=N 



n=N 



< Ojv-1 



l-(l + |t|)e- 



N-y 



N 



\s 



< 



\t\e 



-Ny 



N\ 



\S N \ l-(l + |t|)e- 



-Ny 



In particular, if N > [log(l + 2|t|)]/7, we have < | a jv— 1 1 an d S' N < |ojv-i|/|iSjv|. 
Let us conclude by estimating the size of the largest term max \ a n \. Define 





t 


n | 


v\ 


n(n+l)/2 


(1- 


x\)(l 




\x 


P)-..(l- 


X 


") 



(3.7) 



so that \a n \ < b n (with equality if \t\ = 1 and < x < 1). Suppose that |x| = e~ 7 ; 
it then follows from the computation in ( [3.6|) that 6 n attains its maximum value at 
n = |_l°g(l + |^|)/7j, and that this maximum value is exp[C(|£|)/7 + 0(1)] where 



C(t) = |log(l + t)log 



t 



1 + t 



Li 2 



1 



1 + t 



7T 



(3.8) 



In particular, C(l) = 7r 2 /12 [from ( j2.37[) 1. Therefore, for \t\ = 1 the largest term can 
be as large in magnitude as e 71 " Z 127 (and is indeed of this order when < x < 1); while 
the answer R(t, x) can be as small in magnitude as e~ n Z 67 (and is indeed of this order 
when t = 1 and < x < 1). It is therefore necessary to maintain, in intermediate 
stages of the calculation, approximately (7r 2 /47)/log 10 ~ I.O7/7 digits of working 
precision beyond the number of significant digits desired in the final answer. 



4 Comparison with other algorithms 

Let us conclude by briefly comparing the algorithm based on ( |1.11| ) with some 
alternative algorithms for computing R(t,x). 
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Direct use of the defining product ( |1 . 1| ) manifestly gives an algorithm that is only 
linearly convergent, and in which the convergence rate deteriorates linearly as \x\ j 1. 
Moreover, there is severe loss of numerical precision when multiplying numbers that 
are very near 1. An alternative approach can be based on the logarithmic variant 
Q2.6| ); this sum is again only linearly convergent, but the problem of loss of numerical 
precision is alleviated by use of the logarithm. 

A slight improvement to the algorithm based on ( |1 . 1|) can be obtained by noting 
that 

00 f T N+l 

n = 1 - tz- + o oo - ( 4i ) 



X 

n=N+l 



so that correcting the product ( |1 . 1| ) by the factor 1 — tx N+1 / (1 — x) yields an estimate 
with error 0(x 2N ) rather than 0(x N ). But the basic inefficiencies of the elementary 
algorithm remain. 

Gatteschi |15| has proposed the following iterative algorithm for computing R(t, x):[] 
Choose a complex number a ^ {0,tx} and define 

a = 1 (4.2a) 

A) = — V (4.2b) 
a — tx 

aa n + (1 - o)P n , . . 

a n+1 = a n (4.2c) 

Pn 

Pn+l = «n —T, b- (4-2d) 

xa n + (1 - x)p n 

Gatteschi proves that limn^oo a n = lim^oo (3 n = R(t,x). In fact, it can easily be 
shown by induction that 



a, 



Y[(l - tx k ) (4.3a) 



k=l 



a 



@ n = T^TT a ™ ( 4 - 3b ) 

a — tx n+i 

(though Gatteschi does not note this); so the iteration ( |4.2| ) gives simply a disguised 

way of computing the defining product (|1.1|) and a slight variant of it. Now, it is 
easily seen that 

n f r n+l 

1 + ^— + 0{x 2n ) (4.4a) 



R(t, x) 1 — x 

1 + tx n+1 ( + -) + 0(x 2n ) (4.4b) 



R(t, x) \1 — x a 



9 I have altered his notation to conform to that of the present paper: his a, g, £, x n , y n correspond 
to my tx, x, 1 — <7, a n , Pn- Gatteschi's algorithm has been employed by Allasia and Bonardo 111. 



15 



Therefore, if we set A = 1 + er/(l — x), the linear combination 

atx N+1 

a n = \a n + (1 - \)P n 



(1 - x)(a - tx N + 1 ) 



a, 



(4.5) 



converges to R(t, x) more rapidly than either a n or (3 n does (as Gatteschi observes in 
a special case): namely, a n /R(t,x) = 1 + 0(x 2n ). But this is essentially equivalent 
(modulo higher-order terms) to the "improved" elementary algorithm based on the 
correction factor ( |4.1|) . 

Finally, Slater ^] has computed R(t, x) using the "other" Euler formulaQ 



1 



R(t,x) 



t m x T 



n=l 



m=0 



{l-x){l-x 2 )---{l-x m ) ' 



(4.6) 



But this algorithm is only linearly convergent; it is no better than the logarithmic 
sum (2.6), and indeed is somewhat inferior due to the potentially small denominator. 
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